A nonlinear least squares method for the inverse droplet coagulation problem 
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If the rates, K(x,y), at which particles of size x coalesce with particles of size y is known, then 
the mean-field evolution of the particle-size distribution of an ensemble of irreversibly coalescing 
particles is described by the Smoluchowski equation. We study the corresponding inverse problem 
which aims to determine the coalescence rates, K(x, y) from measurements of the particle size 
distribution. We assume that K(x, y) is a homogeneous function of its arguments, a case which 
occurs commonly in practice. The problem of determining, K(x,y), a function to two variables, 
then reduces to a simpler problem of determining a function of a single variable plus two exponents, 
fj, and v, which characterise the scaling properties of K(x,y). The price of this simplification is 
that the resulting least squares problem is nonlinear in the exponents /i and v. We demonstrate the 
effectiveness of the method on a selection of coalescence problems arising in polymer physics, cloud 
science and astrophysics. The applications include examples in which the particle size distribution is 
stationary owing to the presence of sources and sinks of particles and examples in which the particle 
size distribution is undergoing self-similar relaxation in time. 



I. INTRODUCTION 

Coagulation processes abound in nature and span 
all scales, ranging from the microscopic scales of atmo- 
spheric aerosol formation [l|, to the cosmological scales 
of the clustering of matter within the universe |2| . A par- 
ticular example which played a considerable role in moti- 
vating this work is droplet coalescence in clouds. The role 
of droplet coalescence in the formation and internal dy- 
namics of clouds is of considerable contemporary interest. 
This is because improved understanding of the evolution 
of the droplet size distribution in clouds would increase 
the precision of climate evolution projections Q. Much 
current research in this area focuses on determining the 
rate of coalescence between droplets of different sizes. 
Turbulence in the cloud air mass complicates this task 
significantly. It plays a non-trivial role in determining 
the collision rate of water droplets Direct numer- 

ical simulation of the dynamics of droplets in turbulent 
flows @,H| are possible. It is not clear however that these 
simulations can yet span the range of scales required to 
obtain a full understanding of the role of turbulence in 
facilitating droplet collisions in a cloud Q . For these rea- 
sons an a-priori understanding of droplet collision rates 
in clouds remains elusive. Recent technological advances, 
however, have improved both the quality and quantity 
of empirical data on droplet size distributions at vari- 
ous stages of cloud evolution [icj |. It is therefore timely 
to address the possibility of using these observations to 
solve the inverse problem of determining collision rates 
from measurements of the droplet size distribution. This 
is the topic of this article. We aims to develop a data- 
driven approach to determining collision rates which can 
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complement the insights emerging from current theoret- 
ical and numerical work on this problem. While we have 
motivated our work with in the context of droplet coagu- 
lation in clouds the methods which we develop are quite 
general and provide a quantitative means of constraining 
the choice of model in any coagulation problem in which 
the microphysics is unknown or controversial. 

Throughout this article, we characterise the size of 
droplets by their mass, to, and use N(m, t) to denote 
the droplet size distribution at time t. We denote the 
rate of coalescence between droplets of sizes mi and m,2 
by K(mi,m2)- This function is often referred to as the 
"collision kernel". We assume throughout that the so- 
lution of the forward problem of determining N(m, t) if 
K(m 1,7712) is known is obtained by solving the Smolu- 
chowski coagulation equation (SCE) [ll| : 
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9tN m (t) = 2 / dfn\ K(rn\ , to nii)N mi (t)N m _ mi (t) 
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together subject to an initial condition N m (0) = J\fo(m). 
The final term described a source of particles which in- 
jects "monomers" having mass too at a rate J. Depend- 
ing on the application, J could be zero. We shall take 
the smallest droplet size in the system to be mo = 1. 
We have also explicitly introduced a cutoff mass, M. 
Droplets larger than M are removed from the system. 
Depending on the application, M could be infinite. The 
forward problem is nonlinear in the unknown, N(m,t). 
The inverse problem which forms the topic of this article 
is to determine -FsT(mi,m2) from measurements N(m,t). 
Note that the inverse problem is linear in the unknown, 
-ftf(mi, TO2). The non-triviality of the inverse problem 
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comes from the fact that, since we seek to determine a 
function of two variables from a function of one variable, 
we would generally expect it to be ill-posed. 

Notable previous work on this problem includes the 
work of Onishi and coworkers in the atmospheric science 
context [lH and the work of Ramkrishna and coworkers 
in the chemical engineering context [l3l - [l5j . Onishi et al 
[ijj address the problem of ill-posedness by using signif- 
icant prior knowledge about droplet coalescence in tur- 
bulent conditions to put strong constraints on the func- 
tional form of the kernel. Specifically it was assumed that 
the collision rate could be modeled as a linear superpo- 
sition of the gravitational sedimentations and Saffmann- 
Turner kernels. This simplified the inverse problem to a 
parameter estimation problem at the expense of a loss 
of generality. On the other hand, the methods pio- 
neered by Ramkrishna et al [ll| do not strongly con- 
strain the functional form of K(mi,m2). These au- 
thors address the problem of ill-posedness using a pro- 
cedure known as Tikhonov regularisation. In a previ- 
ous note [HI we explored the ability of the method in 
to solve the inverse problem for kernels of the form, 
K(mi,m2) = \{m\ + wWj) with < A < 1. It was found 
that the method performed relatively poorly when the 
exponent A was fractional. This stems from the fact that 
K{rn\,m2) was represented using Laguerrc polynomials 
which contain only integer powers. 

In this paper, motivated by the fact that many 
practical collision kernels contain fractional powers, we 
present a refined inverse method which deals with the 
problem of fractional exponents up front. Our method 
splits the problem into two stages. In the first stage we 
solve a (nonlinear) parameter estimation problem to de- 
termine a pair of exponents which best describes how 
K(mi,rri2) behaves for large and small masses. In the 
second stage we solve a (linear) inverse problem which 
uses the observations of N(m,t) to correct the detailed 
form of K(nii,m2) without changing the scaling expo- 
nents determined in the first stage. In order to simplify 
our task, we restrict ourselves to cases in which the col- 
lision kernel is a homogeneous function of its arguments. 
While this is common in practice, it is a weakness of our 
approach as compared with that of [l2j which does not 
require this restriction. 

We test our method on two broad classes of inverse 
problems. The first class consists of stationary problems. 
These occur when a source particles is present, J > 
in Eq.((T]), and the sink at M mo is important. For 
sufficiently large times, the droplet size distribution be- 
comes independent of time [l7j and describes a flux of 
mass through the space of droplet sizes from the injec- 
tion scale, mo, to the sink scale, M [l8||. The second 
class consists of time-evolving problems without a source 
of particles in which the droplet size distribution relaxes 
from a prescribed initial condition which we take to be 
monodisperse. For such problems, homogeneous colli- 
sion kernels usually result in the droplet size distribution 
becoming self-similar for large times provided the char- 



acteristic size remains smaller than M . 

The remainder of the article is laid out as follows. 
In Sec|TT]we introduce some properties kinds of the kind 
of collision kernels we expect to find in applications and 
discuss ways of representing such kernels mathematically. 
Next we discuss various aspects of the forward problem 
relevant to our subsequent discussion of the inverse prob- 
lem. The time-dependent forward problem is outlined in 
Scc lIII Al and the stationary forward problem in Sec lIII Bl 
Our main results on the inverse problems are presented in 
Sec lIVI and Sec|V]for the stationary and time-dependent 
inverse problems respectively. Finally, WII presents our 
conclusions and suggestions for further work. 

II. REPRESENTATIONS OF HOMOGENEOUS 
COLLISION KERNELS 

In this paper we focus on scale invariant problems 
for which the kernel, .fT (mi, 7772) is a homogeneous sym- 
metric function of its arguments. We denote the overall 
degree of homogeneity by A: 

K (ami, 0777,2) = a x K(mi, m 2 ). (2) 

Such kernels are important because many physical aggre- 
gation processes exhibit homogeneity for some range of 
scales [1 91 ] - The following model kernel, primarily used in 
the analysis of scaling solutions of the SCE [2(], [2l[ , will 
be of particular importance to us: 

K (m 1 , m 2 ) = I (m^m 2 + m 2 m") . (3) 

where g sets the overall amplitude. Clearly we must have 
A = fx + v. For convenience, we adopt the convention 
throughout that v > /i. The exponents /i and v then 
capture the behaviour of the kernel Eq. (|3]) when one 
mass is much larger than the other: 

Ko(mi, 777,2) ~ m^m^ , with mi <C 777,2. (4) 

Most kernels occuring in practice are not of the form 
([3]). One can however usually [22[ assign a value to the 
exponents /i and v by analogy with Eq. ((3|) by consider- 
ing the behaviour in the limit mi <C 7772. We can then 
write any kernel as a product of (J3j) and another function, 
F(mi,m 2 ): 

K(mi,m 2 ) = Ko(mi,m 2 ) F(mi,m 2 ). (5) 

Since A = /i + v, F must be homogeneous of degree zero. 
It can therefore be expressed as a function, /, of a single 
variable x = mi/m 2 . 

F(mi,m 2 ) = f( 1 ^). (6) 
\m 2 J 

Since K(m\, m 2 ) and, by extension, F(mi, m 2 ) is a sym- 
metric function of its arguments, / must have the sym- 
metry: 

m = f{x- 1 ). (?) 
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Brownian coagulation, Eq.(8) 
Saturn's rings, Eq.(9) 
Differential sedimentation, Eq.(10) 
Nonlinear shear, Eq.(11) 
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FIG. 1. The shape functions, as defined by Eqs.© and ||B). 
for the Brownian coagulation kernel, Eq.([8]), the astrophysical 
coalescence kernel, Eq.©, the differential sedimentation ker- 
nel, Eg. (jlO[) . and the nonlinear shear velocity kernel, Eq.JTT}. 



We will refer to f(x) the shape function of the kernel. 
Since it is a homogeneous function of degree zero it is 
" almost" a constant and by construction must asymptote 
to a constant value as x — > and x — > oo. For the sake of 
concreteness, we have selcted the following kernels from 
the literature to use as test problems in this paper: 
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Eq. ((HJ is the kernel for Brownian coagulation of spherical 
droplets (Tl|. It has v = 1/3 and /x = -1/3. Eq.Q 
is the kernel describing aggregation of ice clusters due 
to differential orbital speed in planetary rings [HI, |2~4| . 
It has v = 2/3 and /i = -1/2. Eq-lfTOj) is the kernel 
most relevant for the cloud problems since it describes 
coalescence of spherical droplets undergoing differential 
sedimentation in the Stokes regime in still air [25[ . It has 
v = 4/3 and fj, = 0. Eq. (fTTj) is the so-called nonlinear 
velocity kernel describing shear-driven coagulation 0, 
[2o| . It has v = 7/9 and // = 0. The respective shape 
functions for each of these kernels are shown in FIGQ] 
Note, how these shape functions are all asymptotically 
constant, and despite the seeming functional complexity 
of the original kernels, have a rather simple form. 

Our approach to solving the inverse problem described 
in the introduction will be to first estimate the expo- 
nents /z and v and then correct the result by an appro- 
priate shape function, f(x). Since the shape function is 



a function of one variable, this should be a considerably 
easier problem. This is very much in the spirit of the 
original work of Ramkrishna and coworkers [IH, [l4[ who 
also exploited the fact that a homogeneous function of 
two variables is determined by its degree and an auxial- 
liary function of a single variable. Our approach is an 
improvement in the sense that the function, f(x), which 
we need to determine is asymptotically constant at large 
and small values making it much easier to deal with. 

To proceed we will need a way of representing func- 
tions on the interval [M _1 , M] (remember we have taken 
mo = 1) which have the symmetry f(x) = f(x~ 1 ). As is 
evident from FigJTJ any symmetric function of log a; has 
the required property. That is we take 



IT l0g(x) 

log(A'/) 
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where h(y) is any symmetric function on the interval 
[— 7T, 7t]. We enforce symmetry indirectly by taking h(y) 
to be a Fourier cosine series truncated after n + 1 terms: 



h(y) 



ao 
2 



fe=i 



a k cos(ky). 



(13) 



One could envisage using other representations but, with 
the possible exception of the differential sedimentation 
kernel, this simple approach will suffice for our purposes. 
Once the exponents [i and v have been determined, the 
inverse problem reduces to using the observed size dis- 
tributions to determine the coefficients of these Fourier 
series. The relative merits of these two approaches will 
be discussed below. 



III. THE FORWARD PROBLEM 

In this section we describe a few features of the so- 
lution of the forward problem which are relevant to our 
subsequent discussion of the inverse problem in the sense 
that we use numerical solutions of Eq.([T]) as input. We 
mention both the time-dependent and stationary cases. 



A. Time-dependent case 

We first consider the evolution of the particle size dis- 
tribution in the absence of a source of monomers. We 
must start from a particular initial distribution, which wc 
usually take to be monodispcrsc: iV m (0) = 5(m — mo). 
If A > 1 then the typical particle size diverges in a fi- 
nite time, t c . This divergence leads to an apparent loss 
of mass from the system as material is absorbed into an 
infinite cluster. This phenomenon is known as gelation 
[27| . Furthermore if v > 1, then gelation occurs instan- 
taneously in the absence of the cut-off, M [28|. In such 
cases, to make sense of Eq.([T]) requires careful considera- 
tion of the regularising role of the cut-off as discussed in 
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FIG. 2. An example of scaling for a time-dependent solution 
of Eq. flj . The main panel shows the time evolution of the par- 
ticle size distribution for the kernel K(rn\,m-2) = v /m i + v /m2. 
The inset panel shows the scaling function, 3>(z), obtained by 
rescaling the data according to Eq. (ll4[) , 



[29(. In order to avoid the complications of gelling sys- 
tems, we shall restrict ourselves here to kernels for which 
A < 1 and v < 1 so that no gelation occurs. 

For homogeneous collision kernels, the time evolution 
of the cluster size distribution tends to become self- 
similar. That is to say, N m (t) tends to the scaling form: 



N(m,t) ~ s(t)- 2 <S>{z) 



s(t)' 



where 



s(t) 



M 2 (t) 
Mi(t) 



(14) 



(15) 



is the typical particle size and "~" denotes the scaling 
limit, m — > oo and s(t) — > oo with z finite. For a modern 
review of the scaling theory of the Smoluchowski equation 
see [2lj|. The scaling function, <&(z) satisfies the following 
equation: 



° = \ J dZl K ( Zlj Z ~ Zl ) $ ( Zl ) $ ( z ~ Zl ) 



$(z) / dzi k(z, zi) $(zi) + \2$(z) + z 



(16) 



Here, k(zi 7 z 2 ) = K(zi,z 2 )/W , where W is the sepa- 
ration constant generated by the self-similarity ansatz 
[20l |2lj . An illustrative example of the kind of data 
which we obtain from a numerical integration of the time- 
dependent forward problem is shown in Fig. [5J 



B. Stationary case 

A stationary cluster size distribution is obtained in the 
limit of large times when a source and sink of particles is 
present. One might imagine that such a stationary state 
could provide a conceptual model of droplet dynamics in 
a cloud where small droplets formed by an ongoing con- 
densation process are driven by air movements to collide 
and coalesce to form larger droplets that eventually be- 
come heavy enough to overcome updrafts and fall out of 
the cloud as rain. The stationary SCE in the presence of 
a source of monomers is: 
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dmi K(mi, m — mi) N rni N r7 



M 



dmi K(m,mi) N mi -\ S(m — mo) (17) 

m 



We distinguish two types of stationary solutions depend- 
ing on whether the kernel has \v — fi\ < 1 or \u — fx\ > 1. 
We refer to the former as "local" kernels and the latter 
nonlocal" for the reasons outlined in UM. For ker- 
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nels having \ v — \i\ < 1 (local), the stationary solution of 
Eq. ([T7j) as M — > 00 takes the power law form 



N„ 



Ay/Jr, 



(18) 



where A is a constant which can be calculated explicitly 
[l7l [l8| . For kernels having |i/ — fi\ > 1 (nonlocal) the 
solution is of the approximate form [30I ] : 



bVjm( t ' 



(19) 



where B is a nonuniversal constant in the sense that it 
depends on M and 7 = v — ^1 — 1. Whether the kernel 
undergoes gelation or not is irrelevant to the stationary 
state. Therefore when we consider the stationary inverse 
problem, we do not need to restrict our choice of kernel to 
the same extent as we do for the time-dependent case. In 
particular we can consider the differential sedimentation 
kernel. 

In order to generate stationary solutions of the SCE we 
can integrate Eq.([T]) forward in time until stationarity is 
achieved. This can be quite slow. Indeed we have shown 
in [3(| that for nonlocal kernels the stationary state can 
become unstable for large M. This curious result means 
that for certain kernels, the stationary state cannot be 
obtained by time integration. For these reasons it is use- 
ful to be able to compute the stationary state directly 
without needing to compute the time evolution numeri- 
cally. The following algorithm achieves this for the model 
kernel Eq.®. The presentation follows the method out- 
lined in the supplementary material of [3p| and is similar 
to the work in [31] . It is included here for completeness. 
To compress the notation it is helpful to introduce the 
moments, 



A4 P , of the size distribution: 



M 



M P =J2 mPN " 



(20) 
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Using the discrete form of equation (|17j) , and a kernel of 
the form in ([3]), we can then use the moments to decom- 
pose (fl~7|) as: 



G + J 



(21) 



where 



m—l 

G = ^2 Ko(mi , m — mi) N mi N, 

mi —1 

2J 



J = —5(m — mo). 
m 



-mi ^22) 

(23) 



Setting mo = 1 gives the stationary monomer density 

2 J 



(24) 



Given the behaviour of the equation for G this permits 
a recursive definition of a stationary distribution, if the 
pair of moments {M.^, M. v ) are known. 

If J, (j, and v are known then (|2ip can be used to in- 
fer the rest of the stationary distribution by treating the 
problem as one of parameter estimation. In this case we 
seek the correct values of the pair of moments (M^, M. v ) 
which will then generate the correct stationary state dis- 
tribution. By treating N m as a function of the pair of 
moments iV m (A^ M , Ai v ) we can create an objective func- 
tion ^{M.^, to be minimised. 

M 

m—l 
M 

+ (M„-J2 m v N m (M^M v )) 2 (25) 

m—l 

(Ma*, M v *) = arg min V(Mu,M v ) (26) 

We remark that this method relies on the special struc- 
ture of the kernel Eq.Q. In general, one cannot avoid 
time integration of the SCE. 



IV. THE STATIONARY INVERSE PROBLEM 

We now present results for the stationary inverse prob- 
lem. The stationary particle size distributions for dif- 
ferent kernels were obtained by numerically solving the 
forward problem either by time integration of Eq.([T]) or 
using the algorithm outlined in Sec. IIII B[ The objective 
is to use the stationary N m to reconstruct K(m,\^m 2 ). 
We take the rate of mass input to be J = 1. In principle, 
the value of J should be obtained from the stationary size 
distribution. For simplicity, we assume that the value of 
J is known a-priori. 



A. Description of the method 

We split the inverse problem into two stages: 

1. Fit the data to the model kernel 

-Ko( TO l, m 2) = 7J (TI1TI2 + ni'2 m l) 

to obtain approximate values for g 7 \x and v. 

2. Estimate the shape function keeping g, /1 and v 
fixed at the values obtained in step 1. 

To implement the first stage, we use the observed val- 
ues of N m to define 

^ m — 1 

T^i(m,g,fi,u) = - ^2 K o{m 1 ,m-m 1 )N mi N m _ mi 

K (m,mi)N mi +6 m ,i. (27) 



mi —1 



M 

£ 

mi — 1 



From these we construct the objective function: 
1 M 

m—l 

We now estimate the values of g, fi and v by minimising 
this function: 

(g*,H*,v*) = arg min S^g, (29) 

(s,M,") 

This is a nonlinear least squares problem since \x and v 
enter into the objective function as exponents. We solve 
it numerically using the Neldcr-Mead algorithm. We 
used the implementation provided in Mathematica™ in 
its Minimize function. 

At the second stage, we fix \i and v to the values ob- 
tained in stage one and introduce a corrected kernel: 

K(m 1 ,m 2 ) = y (mf m 2 + m% m\ ) / {^~^ ( 30 ) 

where the shape function, f(x), is given by Eq. (|T2"]) . It 
depends on n + 1 Fourier coefficients, {afc}£ =0 , through 
Eq. tfTB")) . Choosing an appropriate value for n is impor- 
tant to get good results. If n is too low, Eq. tflUf has 
insufficient flexibility to adequately represent the shape 
function. If n is too high, we start finding implausibly 
oscillatory functions. The choice of n represents our prior 
expectation of how rough or wiggling a function we ex- 
pect the collision kernel to be. In all results presented 
below, the value of n was chosen empirically to be equal 
to 3, giving 4 Fourier coefficients in all. We now proceed 
as before and define 

_^ m—l 

TZ 2 (m, a . . .a n ) = - ^ K{m\,m- m\) N mi N m - mi 

mi— 1 
M 



N m } K(m,mi)N mi + S mA . 



mi — 1 



(31) 
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FIG. 3. Numerical solution of the stationary inverse problem 
for the model kernel, Eq.©, having v = —fj, = 1/3. The 
maximum mass was M = 250 and n + 1 = 4 Fourier coef- 
ficients were used. The numerical values of the parameters 
were v = 0.3333, fx = -0.3332, g = 0.9933, a = 0.9933, 
ai = 0.0268, a 2 = -0.0153 and a 3 = -0.0006. 



From these we construct the objective function: 

M 

M 



1 M 

S 2 (a . ..o n ) = tt X] ^2(m,a . . .o„) 2 



(32) 



The values of the coefficients, arj . . . a n are then obtained 
by solving the linear least squares problem 



(o . 



arg mm 



S 2 (a . . . a n ). 



(33) 



As suggested in [12j , one could weight the objective func- 
tion by the values of N m , on the basis of the heuristic 
that regions with higher density contain more informa- 
tion. We did not find any considerable improvement from 
using such weighting so, in the interest of simplicity, all 
results presented in this paper are unweighted. 



B. Results 

We now show some results obtained by applying this 
method to some stationary size distributions. In order 
to compare the true kernel to the kernels constructed 
from our inverse method we choose to plot slices through 
the kernels as a function of m rather than to show two- 
dimensional surface plots. This is simply to aid clarity. 
We compare a slice through the diagonal, K(m,m), a 
slice through the edge, K(m, 1) and a transverse slice, 
K(M — m,m). Taken together these one-dimensional 
slices give a good overall sense of the quality of the fit. As 
a sanity check, we verified that this algorithm recovers a 
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(a) Brownian coagulation kernel 
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(b) Saturn's rings kernel 
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(c) Nonlinear shear kernel 



250 



FIG. 4. Numerical solutions for the stationary inverse prob- 
lems for the model kernels given in Eqs. (J8j> , ([9| and (|11|) . In 
all cases, the maximum mass was M = 250 and n + 1 = 4 
Fourier coefficients were used 
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reasonable approximation to the kernel if the stationary 
N m is generated using a kernel of the form Eq.© for 
which the shape function is unity. Figure [3] shows the 
results of this test for the case of Eq.© with v = — fi = 
1 /3 and a maximum mass of M = 250. The various slices 
through the kernel mentioned are shown. The points 
representing the reconstructed values closely follow and 
the solid lines representing the true curves. 

We now consider the test kernels, Eq.([51)- (fTT1) which 
have nontrivial shape functions. The values of the expo- 
nents fj, and v at the first stage are usually not partic- 
ularly good. For example, in the case of the Brownian 
coagulation kernel (v = 1/3 and // = —1/3), we obtain 
v* ps 0.24 and /i* ps —0.24. The reason is clear from 
FigfT] with a cut-off of order 10 2 , most of the data are 
in the central region where the shape function is strongly 
varying. Indeed FigfT] suggests that one would need data 
with a cut-off of order 10 6 in order to comfortably en- 
ter the asymptotic regime for of these kernels where the 
"true" values of /i and v should become apparent. 

We then estimated the shape function taking n = 3 
and using Eq. (fT3)) to represent the the shape function. 
The results were not particularly close to the true shape 
function although they seemed qualitatively similar. At 
this point it is important to recognise that neither the 
estimated exponents fx* and v* nor the estimated shape 
function are individually of any importance. What mat- 
ters is the degree to which the combination of the two, as 
written in Eq. (|30p , approximates the true kernel over the 
range of scales for which we have data on the stationary 
size distribution. It is possible for the estimated expo- 
nents fi* and v* and the estimated shape function to be 
individually poor approximations to the "truth" and yet 
provide a good approximation to the true kernel when 
combined. Furthermore, it is ok if the reconstructed ker- 
nel approximates the true kernel very poorly when ex- 
trapolated beyond the range of scales, [1,M], since we 
have no data outside of this range. 

The results for these kernels are shown in Figs. 4(a)| 
|4(b)| and 4(c)| for the Brownian coagulation, Saturn's 
rings and nonlinear shear kernels respectively (see Eqs. 
(JSJ), (j9]) and (jlllO . It is clear that the method does an 
excellent job of recovering an approximation to the truth 
in all cases. The method is therefore quite robust. 



C. The differential sedimentation kernel 

While all seems well at this point, when we applied the 
method to the differential sedimentation kernel, Eq. ([T0"|) . 
most relevant to the cloud physics problem, we found 
that it entirely failed to reconstruct anything reasonable. 
In particular, we obtained negative values for the shape 
function if more than 2 Fourier coefficients were used. We 
can trace the problem to the presence of the cusp in the 
shape function at zero (sec Fig. [TJ). To adequately cap- 
ture this cusp using a Fourier cosine series would require 
the retention of a very large number of terms which, as 



3000 



2500 



2000 



* 1500 



1000 




500 



100 150 200 250 
m 

FIG. 5. Numerical solution for the stationary inverse prob- 
lem for the differential sedimentation kernel, Eq. (|10|l .The 
maximum mass was M = 250 and a fully nonlinear estima- 
tion algorithm with n + 1 = 4 Fourier coefficients was used as 
described in the main text. 



we have already seen, is not a good idea since it provides 
too much freedom to introduce extraneous oscillations. 

In order to get some reasonable results for this kernel, 
after some experimentation, we combined steps 1 and 
2 of the previous method into one single optimization 
problem. This is considerably more expensive numer- 
ically since one now has to do a fully nonlinear mini- 
mization of the objective function in the six-dimensional 
space (v, fj,, oq, oi, 1x3, as). In addition it was necessary 
to include explicit constraints to prevent the shape func- 
tion from becoming negative near 1 which further slow 
down the calculation. The results are shown in Fig. [5] 
While the algorithm recovers the correct qualitative form 
of the kernel, we see there are large quantitative differ- 
ences compared to the results for the other kernels shown 
in Figs 4(a)| |4(b)| and 4(c) Further research will be re- 
quired to improve the performance and reliability of the 
method for such cases. 

Had it not been for the fact that we already knew the 
form of the solution, we probably would not have been 
able to reconstruct the kernel even in the approximate 
way shown in Fig. [5l This reinforces our belief, stated at 
the outset, that the kind of data-driven approach demon- 
strated above can be complementary to existing theoret- 
ical and numerical studies of coagulation phenomena but 
clearly cannot replace them. To conclude our discussion 
of the differential sedimentation kernel, we remark that 
the problems we have encountered due to the cusp in the 
shape function are unlikely to occur if gravitational sed- 
imentation is occuring in a turbulent environment like a 
cloud. This is because spatial variations in the turbulent 
velocity field provide an additional mechanism for colli- 
sions between droplets of equal size as originally pointed 
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FIG. 6. Input data for the time-dependent inverse problem 
in the case of the Brownian coagulation kernel, Eq(S] Main 
panel shows the time evolution of the size distribution and 
the inset shows the data collapse obtained by rescaling this 
data according to Eq. lfHl) . 



to rescale the data according to Eq. (fT4]) . Provided 
that our measured size distributions are in the scal- 
ing regime, this should collapse the data onto a sin- 
gle scaling curve, $(2). This is done for the Brow- 
nian coagulation kernel, Eq.® in Fig. [6] with the 
data collapse shown in the inset. 

2. To discretise Eq. (fT6|) . we need to calculate $(2) 
and its derivative on a regular grid. It is therefore 
convenient to fit the collapsed data to a specific 
functional form. The fit was done using regression 
in the logarithmic variables (u = log.<5>,v = log 2) 
by fitting the the data to the model 

u = bo + biv + b-2V 2 + b^v 3 + b± exp(w) 

and then recovering the required curve by expo- 
nentiation. The result of this fit is superimposed 
on the data in the inset of Figj6l The fitted curve 
can then be differentiated analytically to compute 
the lefthand side of Eq. (fTB|l . 

3. We discretise Eq.fTOl) on Z uniformly spaced 
points, {zi = z m ; n + (i — 1) Az}f =1 where 



out by Saffman and Turner 32[ . This leads to a smooth- 
ing of the cusp and a nonzero value for the collision rates 
of equally sized droplets. In [l2[ it was shown that a 
linear combination of the differential sedimentation and 
Saffman- Turner kernels is a plausible model in the case 
of gravitational settling in turbulent environments. 



V. THE TIME-DEPENDENT INVERSE 
PROBLEM 

In the case of the time-dependent inverse problem, 
some modifications of the method described above per- 
mits the reconstruction of the collision kernel from a suc- 
cession of snapshots of the particle size distribution pro- 
vided that the data span a sufficient range of mass and 
time scales to enter into the scaling regime described in 
section UlI Al The basic idea is to use the scaling ansatz, 
Eq. (fT4|) , to collapse the different snapshots of the parti- 
cle size distribution onto a single scaling function, $(2). 
This curve satisfies Eq. (fT6"|) . Given that $(2) is known, 
the corresponding inverse problem for ^(21,22) is struc- 
turally almost identical to the stationary inverse problem 
which we have already discussed in Sec lIVI After appro- 
priate discretisation, of Eq. p^l) . the methods described 
in Sec. IIVI can be applied with some minor modifications 
of the objective functions Si((j,,v) and ^(ao ■ ■ • to 
take into account the additional linear terms in Eq.(|16p. 

The steps in the procedure are as follows 

I. From the observed snapshots of the size distribu- 
tion, N m (t) we calculate the characteristic parti- 
cle size as given by Eq. ([T5|) and use these values 



(i - \)Az 



and 



A2 = (z max - z min )/(Z - 1). 



2 m in is set by the maximum size reached by s(t). 
We choose 2 max so that it lies in the exponenial 
tail of the scaling function (see inset of FigJ2|) but 
not so large that the $(2) is effectively negligible. 
In the results presented below we typically took 
Z = 250 as the number of discretisation points. 

4. We now proceed as before by defining 

1 4-1 

Ui(i,g,(i,v) = - ^2K (z j ,z l - Zj)$(zj)®(zi - zj) (Az) 
i=i 

z 

- ®( z i) X} Ko(zi> (Az) (34) 

+ 2$(2 l ) + 2 4 — ( Zl ) 

dz 

We again construct the objective function: 

z 



1 Z 



(35) 



and obtain the estimated values of <?, (i and v. 

{g*,V*,v*) = arg min Si{g,n,v). (36) 



9 



5. Finally we correct this result with a shape function: 



(37) 



k(zi,z 2 ) = 9 — (zl z v 2 + z% z\ ) / f ^ 



where the shape function, f(x), is again given by 
Eq. (fT2]) with the n + 1 Fourier coefficients, {afe}fc =0 , 
entering via Eq. (fl~3|) . We now proceed as before and 
define 



i-l 



1 

K 2 (i,a ■ ..a n ) = .- 2J n(zj,Zi - z 3 ) $(zj) ^(z l - Zj) (Az) 

(38) 



3=1 



- $(zi) z j)H z j) (^ z ) 

3 = 1 

+ 2$(zi)+Zi—(zi) 
az 



From these we construct the objective function: 

z 



S 2 (a ...a n ) = — ^ TZ 2 (Mo ■■■««)' 



(39) 



and solving the least squares problem to obtain the 
coefficients: 



(o 



arg mm 

a,Q...a r , 



S 2 (a . ..a n ). 



(40) 



This procedure gives the rescaled kernel, k(z\,z 2 ). To 
return to the original scale requires a knowledge of the 
separation constant, W, which enters when using the 
self-similarity ansatz to split the time-dependent problem 
into an equation for s(t) and a time-independent equa- 
tion for &(z). The presence of this arbitrary constant 
reflects the fact that a rescaling of the amplitude of the 
solution corresponds to a rescaling of time. If we wish 
to fix a value of this constant, we need a way to set the 
time scale. This can be done by fitting the data curve 
obtained using (fTS)) with the analytic solution for s(t): 



»(t) = [(l - x)wt + xy- 



x < 1 



(41) 



Here X is a parameter replacing the initial condition 
term s(0) 1_ \ Provided that (li, v) are retrieved suffi- 
ciently well during the estimation of the kernel, then the 
homogeneity A = fx + v is known and the fitting process 
to obtain W works well. Multiplying the retrieved kernel 
estimate k(zi,z 2 ) by W rescales the result to match the 
original unsealed input kernel K(z\, z 2 ). 

Some illustrative results from the application of this 
method to some time-dependent inverse problems for 
the case of the Brownian coagulation kernel, Eq.©, are 
shown in FigJT] These calculations were done with n = 3. 
Although we get excellent results in this case, in general 
the results tend not to be as good as in the stationary 
case. One reason for this is the increased error made by 
assuming that the system is in the scaling regime. In no 
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FIG. 7. Solution of the time-dependent problem with Brow- 
nian coagulation kernel, Eq.©, based on the data collapse 
shown in Fig[6] The number of discretisation points used was 
Z — 250 and n + 1 = 4 Fourier coefficients were used. The 
numerical results have been rescaled by a factor of 1.8 coming 
from the separation constant (see main text). 



case was the data collapse perfect. A second reason is 
that the exponenially decaying tails of the scaling func- 
tion likely put fewer constraints on the form of the kernel 
since the size distribution becomes negligibly small there. 



VI. CONCLUSIONS AND OUTLOOK 

To conclude, we have presented an approach to solv- 
ing the inverse problem of reconstructing the collision 
kernel from observations of the particle size distribution 
for an ensemble of irreversibly coalescing particles whose 
statistical dynamics is modeled by the Smoluchowski co- 
agulation equation. Compared with previous work on 
this problem, our approach provides a lot of flexibility in 
the functional form of the collision kernel. In particular, 
it handles the possibility of fractional powers of parti- 
cle masses in the kernel, a situation which occurs com- 
monly in applications, in an elegant way. We applied our 
method to a selection of stationary and time-dependent 
inverse problems for which the size distributions were ob- 
tained by numerically solving the forward Smoluchowski 
problem with a variety of different collision kernels taken 
from various branches of physics. The results were of 
sufficiently high quality to demonstrate the feasibility of 
using this kind of data-driven approach to reconstruct 
collision kernels from data. 

Since this problem is under-determined, some prior ex- 
pectation of the kind of functional forms which are rea- 
sonable for the collision kernel is required in order to 
obtain high quality results from our method. This prior 
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expectation is encoded in the choice of the number of 
Fourier coefficients, n, with which to represent the shape 
function of the collision kernel. The choice of n reflects 
the degree- of wiggling which we think is plausible-. It is 
probably possible to come up with cross-validation ar- 
guments to help to select the value of n but it is not 
possible to remove this arbitrariness entirely. For this 
reason, our approach complements rather than replaces 
existing direct theoretical and numerical approaches to 
coagulation phenomena. It is likely to be most useful in 
situations for which quality measurements of the particle 
size distribution are available but for which the underly- 
ing microphysics is unknown or controversial. 

In terms of future research, it would be interesting to 
apply these techniques to some real data. It will there- 
fore be necessary to quantify how well things work in the 
presence of observational noise. It will also be necessary 
to quantify the uncertainty in the collision kernels which 
are reconstructed from the data. To do this, it may be 
best to reformulate the problem in a probabilistic set- 
ting and apply some of the methods of Bayesian inverse 
problems which have been developed recently (see for ex- 
ample [33]). Such a reformulation would also allow us to 



be more explicit about the prior information about the 
kernel which is incorporated into the model and provide 
a means for observational data to over-rule these prior 
choices if they are inconsistent with the observations. 

Finally, the central challenge in the inverse problem is 
that the given information is a function of one variable, 
whereas the kernel is a priori a function of two variables. 
We have shown that constraining K to be homogeneous 
renders the inversion tractable in representative simple 
cases. In doing so we find ourselves using only the steady 
state and/or scaling forms of the cluster size distribution. 
In principle one could exploit the full time dependence 
of the cluster size distribution. This could enable K to 
be inferred as a function of two masses without assuming 
homogeneity. We hope to return to this in future work. 
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